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Abstract. We extend the Kelvin-Helmholtz instability to an expanding background. We study the evolution of a 
non- viscous irrotational fluid and find that for wavelengths much smaller than the Hubble scale small perturbations 
of the fluid are unstable for wavenumbers larger than a critical value. We then apply this result in the early universe, 
treating cold dark matter as a classical fluid with vanishing background pressure. 
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1. Introduction 

In the standard cosmological model quantum fluctuations 
are magnified from microscopic to cosmological scales 
during a period of exponential expansion (inflation) and 
source curvature perturbations with a scale-invariant or 
nearly scale- invariant spectrum ( |Kolb and Turne~1 990 
|Liddle and Lyth, 2000| ). After inflation these perturba- 
tions act as seeds for the large scale structure of the 
universe: first, beginning during radiation domination, 
non-baryonic dark matter "falls" into the potential wells 
formed during inflation and later, during matter domina- 
tion but after decoupling, radiation and baryonic matter 
follow. 

In this picture radiation and baryonic matter can "fall" 
into the potential wells only after decoupling since before 
that their Jeans length is larger than the Hubble radius. 
After decoupling perturbations whose wavelength is larger 
than the Jeans-length grow, but on scales smaller than the 
Jeans-length they oscillate, leaving imprints in the CMB 
that can be observed today. On very small scales structure 
is wiped out through the tight coupling of the baryons 
to the photons until decoupling (Silk, 1968). Unlike in 
baryonic matter and radiation, perturbations in the non- 
baryonic dark matter can start growing as soon as their 
wavelength is smaller than the horizon, since the Jeans 
scale is negligible for non-interacting dark matter. The 
non-baryonic dark matter can therefore cluster and "fall" 
into the potential wells much earlier than the other mat- 
ter. 



In fluid dynamics 
HChandrasekhar, 1981| ) 



is 



the study of instabilities 
well established and has 



also been extended to the astrophysical context (see 
|Inogamov, 19991 and refe rences therein). The Kelvin- 
Helmholtz instability has been known for nearly 150 
years ( |Chand rasckha r7 1981| ). In its simplest form two 
large fluid volumes, separated by an interface, move in 
opposite directions. Neglecting the stabilizing effects of 
density and gravity gradients this setup is unstable for all 
non-zero velocities and wavenumbers. 



We study the initial stages in the evolution of small 
perturbations in an irrotational, non-viscous fluid due to 
a Kelvin-Helmholtz instability in an expanding universe. 
We show that even in such a background these instabilities 
occur. We then apply these results in the early universe. 
Two large colliding regions of cold dark matter might serve 
as an example: the Kelvin-Helmholtz instability can arise 
at the contact surfaces of the two regions. We model the 
dark matter as a fluid with vanishing background pres- 
sure in concordance with the standard cold dark matter 
models | |Kolb and Turner, 1990| |Liddle and Lyth, 2000| > . 
Depending on when this happens in the evolution of the 
background the perturbations arising from the instabil- 
ity grow either exponentially or linearly on scales smaller 
than a critical wavenumber which we determine. On the 
basis of these results we argue that it will be interesting to 
test these perturbations in simulations because they could 
have relevance to the dark matter cusp and halo problems 
and to anisotropies in the CMB. 
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In the next section we give the governing equations for 
the fluid. In Sec. we specify the model geometry and 
the general setup. After giving the background solution in 
Sec. 0J we derive the perturbed equations in Sec. and 
analyse the stability of the system. In Sec. we then cal- 
culate the critical wavenumbers for which instability can 
occur. We apply our results to standard cold dark matter 
in the final section and discuss possible implications. 



Interface C, 




2. Governing equations 

We are interested in the dynamics of a fluid in the early 
universe on scales much smaller than the horizon size, 
i.e. A -C H. We can therefore use a Newtonian ap- 
proximation on an expanding background where only the 
background expansion is described by a flat FRW model 
| |Kolb and Turner, 1990[|Liddle and Lyth, 2000| ). 

We assume that the fluid is non-viscous and irrota- 
tional, its dynamics therefore governed by the Euler equa- 
tion. The fluid does not have to be the only or dominant 
matter component in the model of the universe, we merely 
require that it does not interact with the other fluids such 
as radiation. As outlined above we envisage this fluid to 
describe standard cold dark matter, which has vanishing 
background pressure. We shall nevertheless keep the dis- 
cussion in this and the following sections applicable to the 
general case of a perfect fluid and restrict the analysis 
to cold dark matter only when we consider a particular 
example in the final section. 

The Euler equation in an expanding background is 
given by ^Peebles, 1980| ) 



d 

9? 



Hu l 



d 



1 d 



adx k 



p adx 1 



P 



d 



adx 1 



0, (1) 



where u 1 is the three- velocity of the fluid, P is the pressure, 
p is the energy density of the fluid, <f> is the gravitational 
potential, a — a(t) is the scale factor and H = a/ a is 
the Hubble rate 1 . The energy conservation equation in 
an expanding background is given by 



The condition for the fluid to be irrotational is, 
8 



dx l 



0. 



(2) 



(3) 



3. General formulation 

We can split the variables describing our problem into 
a time-dependent background part and time- and space- 
dependent first order perturbations, 

1> = il>(t) + fy(t,j), (4) 

where ip = p,<f>, P, u % . 

1 Notation: a dot denotes differentiation with respect to co- 
ordinate time t, and i, j — 1,2,3. Comoving length scales A are 
related to the wavenumber by k = 2ir/\. 



Fig. 1. The configuration: fluid moving to the left with 
velocity is separated from fluid moving to the right 
with velocity U a by the perturbed interface 



The geometry of our model is outlined in Fig. ^ two 
large volumes of fluid are separated by an interface at 
x 3 = and move in opposite directions with velocities U a 
and Ub above and below the interface, respectively. We 
allow for a perturbation of the fluid interface separating 
the two flows which is denoted by C = C(*> xl i x2 )- 

The velocity can also be split according to Eq. (T3J, 
u l = U l + u 4 , where XJ % is the background part and v 1 the 
perturbation. Without loss of generality we can choose our 
coordinate system such that the x 1 -coordinate is aligned 
with the direction of the background flow as in Fig. ^ 
The background flow is therefore separated by the a; 3 = 
plane. The background velocity U % has then only one non- 
zero component U 1 = U(t), but for the velocity perturba- 
tion we still have v l — v l (t,x l ). 

We neglect gravitational back reaction on the pertur- 
bations and therefore set 5<j> = 0, i.e. the gravitational 
potential <f> only has a time-dependent zeroth order part. 

4. Background solutions 

For the background described in the previous section we 
get from Eqs. (JIJ and (J2J the simple background Euler 
and energy conservation equations 



U + HU = 0, 
p + 3Hp = , 

which have the solutions 



U = Uo\ — 

a 



P = Po — 



(5) 
(6) 



(7) 



where Uq and po are constants of integration such that at 
t — to we have U = Uq and p = po. 

The background expansion is described by a FRW 
model. We therefore get a power law solution for the scale 
factor ( |Kolb and Turner, 1990||Liddle and Lyth, 2000| ) 



(8) 



where the constant of integration is chosen such that 
a = ciq at t = tg. The exponent 7 characterizes the differ- 
ent epochs of the early and late universe: 7 = 1/2 during 
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radiation domination and 7 = 2/3 during matter domina- 
tion. The evolution of the scale factor is governed by the 
total background energy density. The solution (JSJ) is there- 
fore consistent with our model in the matter dominated 
regime and in the radiation dominated regime despite the 
fact that the dark matter fluid always scales as p oc a~ 3 . 

5. Stability analysis 

In order to analyze the stability of our model we first lin- 
earize the governing equations and then look for growing 
mode solutions. 

Linearizing around the background described in Sec. [31 
we get from Eq. QJ the linearized perturbed Euler equa- 
tion to first order 



■l-u* + Hv* + U 1 - 



d 



1 d 



■5P = 0. 



dt adx 1 p adx 1 

and the energy conservation equation becomes 

1 d 



5p + 3H5p 



a dx L 



^° a dx k 



= 0. 



(9) 



(10) 



The condition for the flow to be irrotational, Eq. J2J), be 
comes 

8 



dx 1 



0. 



(11) 



Substituting Eq. (|f I \ into Eq. I jlOjl , we see that the veloc- 
ity perturbations decouple from the density perturbations. 
To analyse the Kelvin-Helmholtz instability it is therefore 
sufficient to study Eqs. © and lfl"T|) . 

The perturbation in the position of the interface, given 
by x 3 = £(t, x , x 2 ), is a first order quantity and has to be 
included in our analysis. The velocity of the fluid interface 
in the 3-direction is given by 



adx 1 



adx 2 



(12) 



where the RHS is simply the material derivative of the 
physical fluid interface position. Linearizing around the 
background described in Sec. |31 



dt 



HC + U 1 



adx 1 



(13) 



We now make an ansatz for the first order quantities, ex- 
tending Chandrasekhar (Chandrasekhar ,* 1981| ), 



5Q(t,x % ) oc eK^+k^+^t) ( 



(14) 



where ki and k 2 are the components of the wavevector 
and u> is the time dependent frequency. With this ansatz 
growth of the perturbations and hence possible instability 
will set in if ui(t) > 0. Using the ansatz (|I4|I the linearized 
Euler equations JJJJ become 



u> - 



H 



H 



a 

-ill 1 — 



p a 

I d sp 

ap dx 3 



(15) 



where a = 1,2, whereas the continuity equation pif) and 
the linearized kinematic condition Eq. I|13|) give, respec- 
tively, 



d 



i (hv 1 + k 2 v 2 ) + ^ji' 3 = 



u + H + iU 



(16) 



The system of equations l|15JI and Ijlfijl can be solved for 
the velocity perturbation v 3 , 



a J \ dx 3 



(17) 



where k 2 = k\ + k\. We see from Eq. ifTTI) that the x 3 de- 
pendent part of the solution to v 3 has to be oc e ±kx . The 
kinematic condition Ijltjfl shows that v 3 / (u> + H + iU 1 ^) 
has to be continuous across the interface. The solution for 
v 3 is therefore given by 



/•;; = Aiw + H + iU*— ) e' kx \nr.r- , 0. 



for a; 3 < , 



(18) 



where we have used the condition that the velocity has to 
stay finite as x 3 — * ±oo. 

We now integrate Eq. I|17ll across the interface from 
( — e to £ + e and then let e — > to get 



H + ill. 



1 k 1 \ d 

dx 3 



H 



a 



d „ 



(19) 



Combining Eqs. I |18| l and (|19|l finally gives an equation for 
the time derivative of uj, 



uj = —H 

2 a 



i (U & + U h ) ± \U h - U a \ 



(20) 



Since we are studying the growing or decaying behavior of 
the solution we are only interested in the real part of u). 
We can therefore neglect the complex term in the square 
brackets in Eq. H20|) above, because it would only con- 
tribute an oscillatory part to the solution. We also see that 
Eq. I|20[) allows two solutions depending on which sign we 
choose in front of the |Z7b — U a \ term. We only consider 
the "— " sign in the following, since it can already be seen 
that because the Hubble parameter H is always positive 
the "+" choice will only lead to damped solutions. 

Note that we recover the standard result of 
Chandrasekhar {Chandrasekhar, 1981| ), i.e. that all 
wavenumbers are exponentially unstable, if we set a = 1. 

Equation (|20|l can now be integrated using the solution 
for the background velocity Eq. and the scale factor 
Eq. ©. We have \U h - U a \ = U(t/t Q )-~< where we have 
defined U = \ Ubo — Uao | ■ We choose as initial condition 
that there is no amplification at the beginning, i.e. u> = 
at t = to- 
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We study the case of a radiation dominated universe, 
i.e. 7 = 1/2, and a universe with a different exponent 7 
separately. 

Integrating the real part of Eq. <|2U|I we get for the case 
7 = 1/2 



2 V a o 



(21) 



where we have defined t = t/to. 

In the case where 7 ^ 1/2 we get, by integrating the 
real part of Eq. (f^TTf) . 



to = — 7ln£ - 



4(27-1) 
where we have introduced the quantity 
2a 



.4 



(22) 



(23) 



because it plays a crucial role in defining critical wavenum- 
bers in the following section. 



Note that the limit t — > 1 of the RHS of Eq. above is 
well defined, lim^ lnt/(l - i 1 " 27 ) = 1. 

We will now discuss the two cases 7 < 1/2 and 7 > 1 / 2 
separately. 

6.2.1. Case 7 < 1/2 

The case 7 < 1/2 is particularly simple. Using the relation 
( |Abramowitz and Stegun, 1972| ) 



In a; < n{x x l n - 1) for x, n > , 



(26) 



we see that the time dependent part of Eq. Ij25(l is decreas- 
ing and we therefore get, as a conservative value for the 
unstable wavenumber, 



fccrit,L > lA . 



(27) 



For wave numbers satisfying this condition lo > and 
the perturbations grow exponentially with time as ~ 
exp(i 1-27 ) and so the flow is exponentially unstable to 
small perturbations. 



6. Critical wavenumber 

We can now analyze the solutions we have found in the 
previous section and calculate the critical wavenumber be- 
yond which the perturbations are growing. 

6.1. Radiation domination, 7 = 1/2 

In the radiation dominated regime, 7 = 1/2, we find from 
Eq. I|21ll and our ansatz l|14l) that the flow is linearly un- 
stable for all wavenumbers with 



ad 



A 
~2 



(24) 



From Eqs. I|21|l and l|24|) we see that during radiation 
domination perturbations with wavenumbers k\ > fc C rit,r a d 
only grow linearly, i.e, SQ cx t. 

6.2. Case 7 7^ 1/2 

In the case of 7 ^ 1/2 the behavior of the perturbations is 
more complicated than in the radiation dominated case, 
since here the critical wave number is time dependent. 

The solution Eq. (|22|l has two parts with different be- 
haviors: the logarithmic part will always lead to a decrease 
in lo, whereas the second term behaves differently depend- 
ing on whether 7 is smaller or larger than 1/2. In the case 
7 < 1/2 the second term is always increasing and will 
dominate after a finite time for large enough wavenum- 
bers. For 7 > 1/2 the second term in Eq. (12211 will first 
increase and then decrease but can lead to instability if 
it dominates for a long enough period. From Eq. (12211 the 
critical wavenumber for which small perturbations begin 
to grow is 



,(2 7 -l)lnf 
kl > fccrit,7 = lA 



1 - t 1 -^ 



(25) 



6.2.2. Case 7 > 1/2 

From Eq. (|25|l we see that for large i the time dependent 
term lni/(l — t 1 ^ 2 '') tends to Inf. Hence, in this case, 
the critical wave number increases with time. In order to 
get a conservative estimate of the critical wavenumber we 
evaluate Eq. I|25|l at the time when a given perturbation 
is largest. 

The perturbations reach their maximum value at 



"max — tc 



7A 



(28) 



We can now evaluate Eq. I|25ll at t — i max and get an 
implicit equation for the critical wavenumber for 7 > 1/2 



Scrit.G 



In 



vcrit,G 

jA 



(29) 



Solving Eq. I|29|) numerically shows that the error in ne- 
glecting the logarithmic term is ~ 0.1%. We can therefore 
neglect it and get as an estimate that the perturbations 
will grow for wavenumbers 



ki>jA. 



(30) 



When the perturbations reach their maximum they have 
grown as SQ(t mayi ) = <5Q(i )e^ max , where w max = w(i max )- 
Note that we have chosen u>(to) = 0. We therefore find that 
the perturbations have grown at time i max by a factor of 
e Wm ™, where w max is given by 



7 



27 - 1 



jA 



In 



_fci_ 

jA 



1 



(31) 



From Eq. (|26() we see that in Eq. 131|) the square bracket is 
dominated by the first term and hence the perturbations 
are amplified by a factor of ~ e kl ^ A . 
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Our model covers the irrotational, non-viscous case. 
The inclusion of viscosity is even more complicated: for 
small Reynolds numbers it tends to dampen the in- 
stability, whereas for large Reynolds numbers the de- 
velopment of turbulence becomes more likely. Although 
the simplest CDM model does not include viscosity 
there are models, like neutralino dark matter studied in 
(Hof mann et al., 2001| ), which are viscous. We note these 
instabilities will probably not he detected in models based 
on collision-free Boltzmann distributions and that the 
fluid model we use will break down at some small scale 
dependent on the dark matter properties. 

7. Discussion and Conclusions 

In standard fluid dynamics the Kelvin-Helmholtz instabil- 
ity in the incompressible case develops for all wavenum- 
bers k > 0. However in an expanding universe we have 
shown that this is no longer the case and that there are 
critical wave numbers which separate stable and unstable 
domains. These wave numbers depend on the expansion 
of the background and the background velocity. For wave 
numbers larger than the critical wavenumber, k > fc cr it, 
small perturbations will grow. 

More specifically, for a scale factor exponent 7 > 
1/2 perturbations will grow for a wavenumber-dependent 
time. In the case 7 = 1/2, i.e. during radiation domi- 
nation, perturbations grow linearly for all times and for 
7 < 1/2 they grow exponentially for all times. 

As an illustration we look at the case of cold dark mat- 
ter at the beginning of matter domination, 7 = 2/3. The 
physical setting could be two large volumes of cold dark 
matter colliding or some non-homogeneous velocity flows 
in the early universe. In the first case the contact surface 
would be prone to the Kelvin-Helmholtz instability and 
in the second the instability would occur inside the flow. 
From Eq. (|23|) and i|3U|l the physical critical wavenumber 
at the time of matter and radiation equality is given by 

^ = 3L^L ; (32 ) 

where a cq is the scale factor at equality. We choose 
t to be the present age of the universe, t = Hq , 
where H = 100kms~ 1 Mpc~ 1 h, a /a eq = 24000n /t 2 (see 
( |Kolb and Turner, 19901 |Liddle and Lytfa, 2000| )), with 
f2o ~ 1 and the characteristic velocity is U ~ lOOfcms -1 
(see |Green, 2002) |KTyp"in et al., 200 1) . With these num- 
bers the critical scale below which perturbations are am- 
plified is 

A CT it ~ Ikpc. (33) 

Although detailed numerical analysis will be necessary 
in order to make definite quantitative predictions of the 
effect of the instability on the density distribution, we can 
nevertheless speculate that it will lead to a redistribution 
of power from scales ~ X cr it to much smaller scales. 



A significant feature of this example is that the in- 
stability occurs at such a small scale. This is impor- 
tant because it is well known that the cold dark mat- 
ter model has problems at small scales. It seems to 
predict higher densities and more small scale structure 
than is observed. Modifications of the CDM have been 
tried, e.g. Warm Dark Matter and Self-Interacting Dark 
Matter but neither of these solves all the problems, see 
e.g. jPrimack, 2002||Tasitsiomi, 2002| ) for an outline of the 
ideas and detailed references. While these discrepancies 
may be, at least partially, due to numerical and observa- 
tional resolution and selection effects or the consequences 
of baryonic physics, there still appears to be a residual 
problem. 

The surprising clement of the Kelvin-Helmholtz effect 
is that it occurs even on expanding backgrounds and de- 
pending on the value of 7 the instabilities can grow expo- 
nentially. Our results indicate that further analysis using 
numerical methods is needed to explore the effect into the 
nonlinear regime. 
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